'''
Created on Jun 28, 2011

@author: oabalbin
'''

import sys
import numpy as np
import array
import copy
import subprocess
from optparse import OptionParser
from collections import defaultdict,deque #, OrderedDict
from bx.intervals.intersection import Interval, IntervalTree

from exome.vcftools.vcf import *
import exome.vcftools.SNP_positions as snppos 
# For MIONCOSEQ analysis
from exome.fmpn.concensus_snv_creator import main_concesus
#from exome.vcftools.concensus_snv_creator_chron import main_concesus
from exome.vcftools.vcf_bedtools_intersector import main_isec
from exome.variantEval.snps_annotator import create_snps_annotation



if __name__ == '__main__':
    optionparser = OptionParser("usage: %prog [options] ")
    
    optionparser.add_option("--local", dest="local", action="store_true", default=False)
    optionparser.add_option("--cluster", dest="cluster", action="store_true", default=False)
    # Specific options
    # Data files
    optionparser.add_option("-g", "--gatk", dest="gatk",
                            help="gatk vcf file")
    optionparser.add_option("-s", "--sam", dest="sam",
                            help="sam vcf file")
    optionparser.add_option("-v", "--vars", dest="vars",
                            help="vars vcf file")
        
    # Annotation files
    optionparser.add_option("-a", "--actionable", dest="actionable",
                            help="actionable genes bed file")    
    optionparser.add_option("-d", "--dbsnp", dest="dbsnp",
                            help="dbSNP positions file")    
    optionparser.add_option("-p", "--hapmap", dest="hapmap",
                            help="hapmap positions file")
    optionparser.add_option("-m", "--hmgd", dest="hmgd",
                            help="hmgd positions file")
    optionparser.add_option("-i", "--cosmic", dest="cosmic",
                            help="cosmic positions file")


    # Intersector parameters    
    optionparser.add_option("-b", "--bed_tools_path", dest="bed_tools_path",
                            help="bed_tools_path")
    
    optionparser.add_option("-l", "--ref_probe_file",dest="ref_probe_file",
                            help="reference file with probes or exons")
    
    optionparser.add_option("-t", "--vcf_file_annotated",dest="vcf_file_annotated", 
                            action="store_true", default=False,
                            help="Boolean that indicates vcf_file is annotated")

        
    # Final Output Report
    optionparser.add_option("-o", "--output_file", dest="output_file",
                            help="output_file")
    optionparser.add_option("-c", "--complement", action='store_true', default=False,
                            help="output_file")

    
    (options, args) = optionparser.parse_args()
    print options
    
    # Fixed variables to tun in exds machine:
    path_to_vcfCodSNPs = '/exds/users/oabalbin/sw/vcfCodingSnps.v1.5/'
    reference_genome='/exds/users/oabalbin/sw/vcfCodingSnps.v1.5/referenceGenomes/human_g1k_v37.fasta'
    genelist_annotation='/exds/users/oabalbin/sw/vcfCodingSnps.v1.5/geneLists/refGene.B37.txt'
    path_to_exons_bed_file='/exds/projects/alignment_indexes/exome_targets/agilent/SureSelect_50Mb_Kitv3/refGene.B37.exons.pad50.bed'
    path_to_bedtools= '/exds/users/oabalbin/sw/bedtools/current/bin/'
        
    # Input files:
    vcf_files, annot_files,isec_files_dict,annot_vcf_files = \
    defaultdict(),defaultdict(),defaultdict(),defaultdict()
    
    #if not (options.gatk ^ options.sam ^ options.vars):
    #    print >> sys.stderr, "You need a minimum of 1 vcf file to do the analysis"
    #    sys.exit(0)
    if options.gatk:
        vcf_files['GATK']=options.gatk
    if options.sam:
        vcf_files['SAM']=options.sam
    if options.vars:
        vcf_files['VARS']=options.vars
    if options.vars=="" and options.gatk=="" and options.sam=="":
        print >> sys.stderr, "You need a minimum of 1 vcf file to do the analysis"
        sys.exit(0)

        
    if options.dbsnp:
        annot_files["dbSNP"]=options.dbsnp
    if options.hapmap:
        annot_files["HapMap"]=options.hapmap
    if options.hmgd:
        annot_files["HMGD"]=options.hmgd
    if options.cosmic:
        annot_files["COSMIC"]=options.cosmic
    if options.actionable:
        annot_files["Actionable"]=options.actionable
    if options.dbsnp=="" and options.hapmap=="" and options.hmgd=="" and options.actionable=="":
        print >> sys.stderr, "You need a minimum of 1 annotation file to do the analysis"

    
    
    
    # Annotate the vcf if it is not annotated
    
    if not options.vcf_file_annotated:
        for name, input_vcf_file in vcf_files.iteritems():
            output_vcf_file=input_vcf_file.replace('.vcf','_annot.vcf')
            output_vcf_log = input_vcf_file.replace('.vcf','_annot.log')
            annot_vcf_files[name]=output_vcf_file
            # This can be in an external function
            subprocess.call(create_snps_annotation(input_vcf_file, output_vcf_file, output_vcf_log,
                                   reference_genome, genelist_annotation, path_to_vcfCodSNPs))
    else:
        annot_vcf_files=copy.copy(vcf_files)
    
    
    # Filtering files to exon or pull down regions only.
    local=True
    if not options.bed_tools_path:
        bedtools_path=path_to_bedtools
    else:
        bedtools_path=options.bed_tools_path
            
    if not options.ref_probe_file:
        fileB=path_to_exons_bed_file
    else:
        fileB=options.ref_probe_file

    # Do the intersections using bedtools.
    for name, query_vcf_file in annot_vcf_files.iteritems():
        isecfile = query_vcf_file.replace('.vcf','_isec_exomePad50.vcf')
        hederfile = query_vcf_file.replace('.vcf','_header.txt')
        isecfile_header=isecfile.replace('.vcf','h.vcf')
        
        isec_files_dict[name]=isecfile_header
        
        main_isec(bedtools_path, query_vcf_file, fileB, isecfile, 
                  isecfile_header, hederfile,local)

    
    
    # Create Concensus Report

    outfile =options.output_file+'.txt'
    mmapfile =options.output_file+'.memmap'
    print options
    
    main_concesus(isec_files_dict, annot_files, outfile, mmapfile, options.complement)
